module soil_carbon_mod

!Don't use external model stuff if compiling standalone version
#ifndef STANDALONE_SOIL_CARBON

#include "../shared/debug.inc"

use constants_mod, only : pi, dens_h2o
use land_constants_mod, only : Rugas, seconds_per_year
use fms_mod, only: check_nml_error, file_exist, close_file, &
            stdlog, mpp_pe, mpp_root_pe, error_mesg, FATAL, NOTE
use vegn_data_mod, only: K1,K2
use land_numerics_mod,only: tridiag
use land_debug_mod, only: get_current_point, is_watch_point
use land_data_mod, only: log_version

#ifdef INTERNAL_FILE_NML
use mpp_mod, only: input_nml_file
#else
use fms_mod, only: open_namelist_file
#endif
#endif


!End of external model stuff

implicit none

private


! ==== public interfaces =====================================================
public :: soil_carbon_pool
public :: soilMaxCohorts

public :: update_pool
public :: add_litter
public :: add_carbon_to_cohorts
public :: combine_pools
public :: poolTotalCarbon
public :: pool_total_carbon
public :: init_soil_carbon
public :: read_soil_carbon_namelist
public :: remove_carbon_fraction_from_pool
public :: litter_leaching
!public :: carbon_leaching
public :: carbon_leaching_with_litter
public :: cull_cohorts
public :: transfer_pool_fraction
public :: retrieve_DOC ! report DOC concentration to hlsp_hydrology
#ifndef STANDALONE_SOIL_CARBON
public :: adjust_pool_ncohorts
public :: A_function
#endif
public :: debug_pool

! public data
public :: soil_carbon_option
! =====end of public interfaces ==============================================


! ==== module constants ======================================================
integer, parameter, public :: N_C_TYPES=3  !Carbon chemical species (Cellulose, lignin, microbial products)
integer, parameter, public :: & ! indices of carbon chemical species
    C_CEL = 1, & ! cellulose (fast)
    C_LIG = 2, & ! lignin (slow)
    C_MIC = 3    ! microbial producs

character(len=*), parameter :: module_name = 'soil_carbon_mod'
#include "../shared/version_variable.inc"

integer, parameter, public :: &
    SOILC_CENTURY          = 1, &
    SOILC_CENTURY_BY_LAYER = 2, &
    SOILC_CORPSE           = 3
integer, parameter :: init_n_cohorts = 3 ! initial number of cohorts in a litter pool

#ifdef STANDALONE_SOIL_CARBON
real,parameter::pi=3.141592
real,parameter::Rugas=8.314472
real,parameter::dens_h2o=1000.
integer,parameter::FATAL=0,NOTE=1
real,parameter :: seconds_per_year=86400.0*365.0
#endif

! ==== types =================================================================

!Individual litter cohort with its own carbon pools
type litterCohort
    !Carbon pools
    real,dimension(n_c_types)::litterC  !Carbon substrate. Relative amounts of C species are important for resp rate
    real,dimension(n_c_types)::protectedC !Aggregate and mineral complex C
    real::livingMicrobeC               !Current carbon mass of live microbes
    real::CO2                          !Cumulative CO2 generated by decomposition
    real::Rtot                         !Cumulative decomposition (includes double counting through microbial products)
    !Original C
    real::originalLitterC              !Keep track for carbon balance check
end type litterCohort


!Pool type, can hold multiple individual litter/carbon cohorts
type soil_carbon_pool
    integer::max_cohorts    !Maximum number of cohorts, allows different pools to have different maximum heterogeneity
    integer::n_cohorts = 0  !actual number of cohorts in the pool
    real::protection_rate   !Pool-specific rate that carbon is transferred to protected pool (aggregate and chemical protection)
    real::Qmax              !Pool DOC sorption capacity (See Mayes et al 2012)
    type(litterCohort),allocatable::litterCohorts(:)
    real::dissolved_carbon(n_c_types)
end type soil_carbon_pool

!==== module variables =======================================================

!---- namelist ---------------------------------------------------------------
character(32) :: soil_carbon_model_to_use = 'CENTURY-like' ! or 'CENTURY-like-by-layer', or 'CORPSE'
real,dimension(n_c_types) :: Ea=(/37e3,54e3,50e3/)          !Activation energy (kJ/mol)
real,dimension(n_c_types) :: vmaxref=(/4e1,1e1,.5e1/)       !Vmax at reference temperature (yr-1)
real,dimension(n_c_types) :: kC=(/.5,.1,0.05/)              !Michaelis-Menton C parameter (dimensionless)
real                      :: Tmic=0.2                       !Microbial turnover rate (yr-1)
real                      :: et=0.5                         !Fraction of microbial turnover not converted to CO2
real,dimension(n_c_types) :: eup=(/0.6,0.2,0.1/)            !Fraction of respired C that goes into microbial biomass
real                      :: minMicrobeC=1e-5               !Minimum microbial biomass (prevents complete collapse if > 0.0)
real                      :: combineThreshold=1e-2          !Controls when cohorts are combined
real                      :: D=3.0                          !Diffusion coefficient (for carbon and enzyme movement)
real                      :: sol=0.4                        !Substrate solubility (for carbon access by microbes)
real                      :: enzfrac=1.0                    !Relative amount of enzymes produced by microbes
real                      :: tProtected=10.0                !Turnover rate of protected carbon (yr-1)
real                      :: protection_rate=1.0            !Rate that carbon becomes protected (yr-1 kg-microbial-biomass-1)
real,dimension(n_c_types) :: protection_species=(/0.5,0.5,1.0/)  !Relative protection rate of each flavor
real                      :: leaching_solubility=1e-2       !Rate carbon dissolves in soil water at saturated moisture (yr-1)
real,dimension(n_c_types) :: flavor_relative_solubility=(/1.0,1.0,1.0/) !For each C flavor, relative to 1.0
real                      :: protected_relative_solubility=1.0 !Relative to 1.0
real                      :: DOC_deposition_rate=1e-2       !Rate carbon is deposited from DOC ((kgC/kgH2O)-1 yr-1)
real                      :: gas_diffusion_exp=2.5          !Exponent for gas diffusion power law dependence on theta
                                                            !See Meslin et al 2010, SSAJ
real                      :: litterDensity=22.0             !C density of litter layer (kg/m3)
                                                            !22.0 roughly from Gaudinsky et al 2000
real                      :: min_anaerobic_resp_factor=0.0  !Minimum for high soil moisture Resp limitation

integer                   :: soilMaxCohorts=15              !Maximum number of cohorts in soil carbon pools
real                      :: tol=1e-4                       !Tolerance for cohort carbon check

namelist /soil_carbon_nml/ &
            soil_carbon_model_to_use, &
            Ea,vmaxref,kC,Tmic,et,eup,minMicrobeC,combineThreshold,soilMaxCohorts,gas_diffusion_exp,&
            tol,enzfrac,tProtected,protection_rate,protection_species,leaching_solubility,flavor_relative_solubility,DOC_deposition_rate,&
            litterDensity,protected_relative_solubility,min_anaerobic_resp_factor

!---- end-of-namelist --------------------------------------------------------
integer, protected :: soil_carbon_option = 0    ! flag specifying which soil carbon to use,
        ! one of SOILC_CENTURY, SOILC_CENTURY_BY_LAYER, SOILC_CORPSE


contains ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-


subroutine init_soil_carbon(pool,protectionRate,Qmax,max_cohorts)
    type(soil_carbon_pool),intent(inout)::pool
    real,optional,intent(in) :: protectionRate,Qmax
    integer,optional,intent(in) :: max_cohorts

    pool%max_cohorts=soilMaxCohorts
    pool%protection_rate=protection_rate
    pool%Qmax=0.0
    pool%dissolved_carbon=0.0

    IF (present(max_cohorts)) pool%max_cohorts=min(max_cohorts,soilMaxCohorts)
    IF (present(protectionRate)) pool%protection_rate=protectionRate
    IF (present(Qmax)) pool%Qmax=Qmax

    pool%n_cohorts = 0
    IF (allocated(pool%litterCohorts)) deallocate(pool%litterCohorts)
end subroutine init_soil_carbon

! =============================================================================
#ifndef STANDALONE_SOIL_CARBON
subroutine read_soil_carbon_namelist
  integer :: unit         ! unit for namelist i/o
  integer :: io           ! i/o status for the namelist
  integer :: ierr         ! error code, returned by i/o routines
  integer :: i
  real    :: z

  call log_version(version, module_name, &
  __FILE__)
#ifdef INTERNAL_FILE_NML
  read (input_nml_file, nml=soil_carbon_nml, iostat=io)
  ierr = check_nml_error(io, 'soil_carbon_nml')
#else
  if (file_exist('input.nml')) then
     unit = open_namelist_file()
     ierr = 1;
     do while (ierr /= 0)
        read (unit, nml=soil_carbon_nml, iostat=io, end=10)
        ierr = check_nml_error (io, 'soil_carbon_nml')
     enddo
10   continue
     call close_file (unit)
  endif
#endif
  if (mpp_pe() == mpp_root_pe()) then
     unit=stdlog()
     write(unit, nml=soil_carbon_nml)
  endif


  ! parse soil carbon option
  select case (soil_carbon_model_to_use)
  case('CENTURY-like')
    soil_carbon_option = SOILC_CENTURY
  case('CENTURY-like-by-layer')
    soil_carbon_option = SOILC_CENTURY_BY_LAYER
  case('CORPSE')
    soil_carbon_option = SOILC_CORPSE
  case default
    call error_mesg('read_soil_carbon_namelist', &
        '"'//trim(soil_carbon_model_to_use)//'" is an invalid option for soil_carbon_model_to_use', FATAL)
  end select

end subroutine read_soil_carbon_namelist
#endif


#ifdef STANDALONE_SOIL_CARBON
subroutine read_soil_carbon_namelist(file)
character*(*),optional,intent(in)::file

integer :: namelistunit
namelistunit=12300
!Read parameters from namelist
if(present(file)) then
print *,'Reading soil params from namelist: ',file
OPEN(unit=namelistunit,file=file)
else
OPEN(unit=namelistunit,file='soilparams.nml')
endif
READ(unit=namelistunit,NML=soil_carbon_nml)
CLOSE(unit=namelistunit)

end subroutine
#endif


subroutine update_pool(pool,T,theta,air_filled_porosity,liquid_water,frozen_water,dt,layerThickness,&
            fast_C_loss_rate,slow_C_loss_rate,deadmic_C_loss_rate,CO2prod,deadmic_produced,protected_produced,protected_turnover_rate,&
            C_dissolved,deposited_C,badCohort)
    type(soil_carbon_pool),intent(inout)::pool
    real,intent(in)::T,theta,dt,air_filled_porosity,liquid_water,frozen_water,layerThickness
    real,intent(out)::fast_C_loss_rate,slow_C_loss_rate,deadmic_C_loss_rate,CO2prod,&
            protected_produced(n_c_types),protected_turnover_rate(n_c_types),deadmic_produced
    real,intent(out)::C_dissolved(n_c_types),deposited_C(n_c_types)
    integer,intent(out),optional::badCohort
    ! dt is in years!

    integer::n
    real::tempresp(n_c_types),resp(n_c_types),temp_protected(n_c_types),&
            tempCO2,poolProtectedC,temp_protected_turnover_rate(n_c_types),&
            Prate_limited(n_c_types),prevC(n_c_types),temp_deadmic
    real::dissolution_rate,liquid_frac
    real::protected_dissolved(n_c_types),livemicrobe_dissolved
    real::cohortVolume(pool%max_cohorts),activeVolume,inactiveVolume
    real::protected_solubility

!   if (is_watch_point()) then
!      write(*,*)'##### update_pool input ####'
!      __DEBUG4__(T,theta,dt,air_filled_porosity)
!      __DEBUG3__(liquid_water,frozen_water,layerThickness)
!      do n=1, pool%n_cohorts
!          write(*,*) 'cohort ',n
!          __DEBUG1__(pool%litterCohorts(n)%litterC)
!          __DEBUG1__(pool%litterCohorts(n)%protectedC)
!           __DEBUG2__(pool%litterCohorts(n)%livingMicrobeC,pool%litterCohorts(n)%originalLitterC)
!           __DEBUG1__(pool%litterCohorts(n)%CO2)
!       enddo
!    endif

    resp=0.0
    CO2prod=0.0
    protected_produced=0.0
    protected_turnover_rate=0.0
    deadmic_produced=0.0

    if(present(badCohort)) badCohort=0

    if(.NOT.allocated(pool%litterCohorts)) call add_litter(pool,(/0.0,0.0,0.0/))
    call cull_cohorts(pool)

    !Update dissolved carbon
    !If there is liquid water, it exchanges carbon with the soil
    !If all the water is frozen, there is no exchange
    !If there is no water at all, then any dissolved carbon is dumped into the soil
    if(liquid_water > 0.0) then
        liquid_frac=liquid_water/(liquid_water+frozen_water)
        deposited_C=DOC_deposition_rate*pool%dissolved_carbon*liquid_frac*dt
        dissolution_rate=leaching_solubility*theta
        where(deposited_C > pool%dissolved_carbon)
            deposited_C=pool%dissolved_carbon
            pool%dissolved_carbon=0.0
        elsewhere
            pool%dissolved_carbon=pool%dissolved_carbon-deposited_C
        end where

        !Protected carbon can dissolve, but much faster under high moisture conditions
        protected_solubility=theta**gas_diffusion_exp*protected_relative_solubility
        if (protected_solubility<0.0) protected_solubility=0.0
        if (protected_solubility>1.0) protected_solubility=1.0

        call remove_carbon_fraction_from_pool(pool,dissolution_rate*dt,&
                litter_removed=C_dissolved,protected_removed=protected_dissolved,&
                livemicrobe_removed=livemicrobe_dissolved,&
                protectedMobility=protected_solubility,livingMicrobeMobility=0.0,&
                litterMobility=flavor_relative_solubility)

        C_dissolved=C_dissolved+protected_dissolved
        protected_turnover_rate=protected_turnover_rate+protected_dissolved/dt
    elseif(frozen_water > 0.0) then
        deposited_C=0.0
        C_dissolved=0.0
    else
        deposited_C=DOC_deposition_rate*pool%dissolved_carbon*dt
        where(deposited_C > pool%dissolved_carbon)
            deposited_C=pool%dissolved_carbon
            pool%dissolved_carbon=0.0
        elsewhere
            pool%dissolved_carbon=pool%dissolved_carbon-deposited_C
        end where
        C_dissolved=0.0
    endif

    call add_carbon_to_cohorts(pool,litterC=deposited_C)
    pool%dissolved_carbon=pool%dissolved_carbon+C_dissolved



    !Protection rate is multiplied by available space, so protected C does not exceed Pmax
    !However, the rate is not adjusted for each cohort, so it may go slightly above the pool Pmax
    ! ---   Based on conversation with Melanie Mayes, I am changing this so Qmax affects the protected
    ! ---   carbon formation rate rather than the maximum.
!    call poolTotalCarbon(pool,protectedC=poolProtectedC)
!    IF (pool%Pmax.eq.0.0) THEN
!        Prate_limited=0.0
!    ELSE
!        Prate_limited=max(0.0,pool%protection_rate*protection_species*(1.0-poolProtectedC/pool%Pmax))
!    ENDIF
    Prate_limited=pool%protection_rate*protection_species*pool%Qmax

    !Need to convert originalLitterC into a meaningful volume, since it will keep increasing as cohorts are combined
    !How about this: non-mineralized C keeps volume based on estimated density
    !Volume of mineralized C is just capped at remaining layer volume, with the assumption that the mineralized portion
    !of volume for all cohorts just gets intermingled and doesn't need to sum to layer volume
!    if (is_watch_point()) then
!       write(*,*) '##### update_pool outpt #####'
!    endif
    do n=1, pool%n_cohorts
!        if (is_watch_point()) then
!           write(*,*) 'cohort ',n
!           __DEBUG1__(pool%litterCohorts(n)%litterC)
!           __DEBUG1__(pool%litterCohorts(n)%protectedC)
!           __DEBUG2__(pool%litterCohorts(n)%livingMicrobeC,pool%litterCohorts(n)%originalLitterC)
!           __DEBUG1__(pool%litterCohorts(n)%CO2)
!        endif
        activeVolume=cohortCsum(pool%litterCohorts(n),.TRUE.)/litterDensity
        inactiveVolume=min(pool%litterCohorts(n)%originalLitterC/litterDensity,layerThickness)-activeVolume
        cohortVolume(n)=activeVolume+max(0.0,inactiveVolume)
    enddo


    DO n=1,pool%n_cohorts
        prevC=pool%litterCohorts(n)%litterC
        call update_cohort(pool%litterCohorts(n),cohortVolume(n),T,max(theta,0.0),max(air_filled_porosity,0.0),&
                        Prate_limited,dt,tempresp,temp_deadmic,temp_protected,temp_protected_turnover_rate,tempCO2)
        IF (.NOT. check_cohort(pool%litterCohorts(n))) THEN
            if(present(badCohort)) badCohort=n
            WRITE (*,*) 'UPDATE_POOL: Cohort',n,'of',pool%n_cohorts,'bad'
            call print_cohort(pool%litterCohorts(n))
            WRITE (*,*) 'Dissolved carbon =',pool%dissolved_carbon
            WRITE (*,*) 'Latest respiration:',tempResp*dt
            WRITE (*,*) 'Previous unprotected C:',prevC
        ENDIF

        resp=resp+tempresp
        protected_produced=protected_produced+temp_protected
        protected_turnover_rate=protected_turnover_rate+temp_protected_turnover_rate
        deadmic_produced=deadmic_produced+temp_deadmic
        CO2prod=CO2prod+tempCO2
!        if (is_watch_point()) then
!           __DEBUG4__(cohortVolume(n),T,theta,air_filled_porosity)
!           __DEBUG3__(Prate_limited,tempCO2,CO2prod)
!        endif
    ENDDO

    fast_C_loss_rate=resp(1)
    slow_C_loss_rate=resp(2)
    deadmic_C_loss_rate=resp(3)
end subroutine update_pool


!Do litter respiration and microbial turnover for one litter cohort
subroutine update_cohort(cohort,cohortVolume,T,theta,air_filled_porosity,protection_rate,dt,totalResp,deadmic_produced,protected_produced,protected_turnover_rate,CO2prod)
    type(litterCohort),intent(inout)::cohort
    real,intent(out)::totalResp(n_c_types),CO2prod,protected_produced(n_c_types),protected_turnover_rate(n_c_types),deadmic_produced
    real,intent(in)::cohortVolume,T,theta,dt,protection_rate(n_c_types) !Temperature, volumetric water content, delta time
    real,intent(in)::air_filled_porosity

    real::microbeTurnover,temp_microbes
    real::tempResp(n_c_types),oldminC(n_c_types)
    real,dimension(n_c_types)::protectedCTurnover,newProtectedC

    totalResp=0.0
    tempResp=0.0
    CO2prod=0.0
    deadmic_produced=0.0

    !Decompose carbon pools
    !Litter
    tempResp=Resp(cohort%litterC,cohort%livingMicrobeC,T,theta,air_filled_porosity)
    where(dt*tempResp > cohort%litterC)
        tempResp=cohort%litterC/dt
        cohort%litterC=0.0
    elsewhere
        cohort%litterC=cohort%litterC-dt*tempResp
    end where


    totalResp=totalResp+tempResp


    !Update microbial biomass
    temp_microbes=cohort%livingMicrobeC
    microbeTurnover=max(0.0,(cohort%livingMicrobeC-minMicrobeC*sum(cohort%litterC))/Tmic)
    cohort%livingMicrobeC=cohort%livingMicrobeC + dt*(sum(eup*totalResp) - microbeTurnover)
    deadmic_produced=dt*microbeTurnover*et
    cohort%litterC(3)=cohort%litterC(3)+deadmic_produced
    CO2prod=dt*(sum(totalResp*(1.0-eup))+microbeTurnover*(1.0-et))
    IF(cohort%livingMicrobeC<0) then
        print *,'Live microbe C',cohort%livingMicrobeC
        print *,'Respired C',totalResp*dt
        print *,'microbe turnover',dt*microbeTurnover
        print *,'Previous microbe C',temp_microbes
        print *,'dt',dt
        call print_cohort(cohort)
        call error_mesg('update_cohort','Microbe C < 0',FATAL)
    ENDIF
    cohort%CO2=cohort%CO2+CO2prod


    !Update protected carbon
    protectedCturnover = cohort%protectedC/tProtected
    IF(sum(cohort%litterC).gt.0.0 .and. cohortVolume.gt.0.0) THEN
        !Change: divide by volume instead of litter C. Keeps it linear with size, but allows dependence on unprotected C
        newProtectedC = protection_rate*cohort%livingMicrobeC/cohortVolume*cohort%litterC*dt
    ELSE
        newProtectedC = 0.0
    ENDIF
    where(newProtectedC.gt.cohort%litterC) newProtectedC=cohort%litterC
    cohort%protectedC = cohort%protectedC + newProtectedC - dt*protectedCturnover
    cohort%litterC = cohort%litterC - newProtectedC + dt*protectedCturnover

    protected_produced=newProtectedC
    protected_turnover_rate=protectedCturnover

    cohort%Rtot=cohort%Rtot+sum(totalResp)
end subroutine update_cohort





! COPIED FROM VEGN_DYNAMICS
! ============================================================================
! The combined reduction in decomposition rate as a funciton of TEMP and MOIST
! Based on CENTURY Parton et al 1993 GBC 7(4):785-809 and Bolker's copy of
! CENTURY code
elemental function A_function(soilt, theta) result(A)
  real :: A                 ! return value, resulting reduction in decomposition rate
  real, intent(in) :: soilt ! effective temperature for soil carbon decomposition
  real, intent(in) :: theta

  real :: soil_temp; ! temperature of the soil, deg C
  real :: Td; ! rate multiplier due to temp
  real :: Wd; ! rate reduction due to mositure

  ! coefficeints and terms used in temperaturex term
  real :: Topt,Tmax,t1,t2,tshl,tshr;

  soil_temp = soilt-273.16;

  ! EFFECT OF TEMPERATURE
  ! from Bolker's century code
  Tmax=45.0;
  if (soil_temp > Tmax) soil_temp = Tmax;
  Topt=35.0;
  tshr=0.2; tshl=2.63;
  t1=(Tmax-soil_temp)/(Tmax-Topt);
  t2=exp((tshr/tshl)*(1.-t1**tshl));
  Td=t1**tshr*t2;

  if (soil_temp > -10) Td=Td+0.05;
  if (Td > 1.) Td=1.;

  ! EFFECT OF MOISTURE
  ! Linn and Doran, 1984, Soil Sci. Amer. J. 48:1267-1272
  ! This differs from the Century Wd
  ! was modified by slm/ens based on the figures from the above paper
  !     (not the reported function)

  if(theta <= 0.3) then
     Wd = 0.2;
  else if(theta <= 0.6) then
     Wd = 0.2+0.8*(theta-0.3)/0.3;
  else
     Wd = exp(2.3*(0.6-theta));
  endif

  A = (Td*Wd); ! the combined (multiplicative) effect of temp and water
               ! on decomposition rates
end function A_function


subroutine initializeCohort(cohort,litterInputC,initialMicrobeC,initialProtectedC,multiplier,CO2)
    type(litterCohort),intent(out)::cohort
    real,intent(in),optional::litterInputC(n_c_types),CO2
    real,intent(in),optional::multiplier
    real,intent(in),optional::initialMicrobeC
    real,intent(in),optional::initialProtectedC(n_c_types)

    real::litterInputCval(n_c_types),initialMicrobeCval,CO2val,multiplierval,initialProtectedCval(n_c_types)

    litterInputCval=0.0
    !initialMinC=0.0
    initialProtectedCval=0.0
    initialMicrobeCval=0.0
    CO2val=0.0
    multiplierval=1.0

    if (present(litterInputC)) litterInputCval=litterInputC
    !if (present(initialMinC)) initialMinCval=initialMinC
    if (present(initialprotectedC)) initialProtectedCval=initialProtectedC
    if (present(initialMicrobeC)) initialMicrobeCval=initialMicrobeC
    if (present(CO2)) CO2val=CO2
    if (present(multiplier)) multiplierval=multiplier


    !originalLitterC is there to check carbon conservation
    cohort%originalLitterC=(sum(litterInputCval)+CO2val+initialMicrobeCval+sum(initialProtectedCval))*multiplierval
    cohort%litterC=litterInputCval*multiplierval
    cohort%protectedC=initialProtectedCval*multiplierval
    !cohort%minC=initialMinCval*multiplierval
    cohort%livingMicrobeC=initialMicrobeCval*multiplierval
    cohort%CO2=CO2val*multiplierval
    cohort%Rtot=0.0

end subroutine initializeCohort


function Resp(Ctotal,Chet,T,theta,air_filled_porosity)
    real,intent(in)::Chet                       !heterotrophic (microbial) C
    real,intent(in)::T,theta                    !temperature (k), theta (fraction of 1.0)
    real,intent(in)::air_filled_porosity        !Fraction of 1.0.  Different from theta since it includes ice
    real,intent(in),dimension(n_c_types)::Ctotal !Substrate C
    real,dimension(n_c_types)::Resp,tempresp
    real::enz,Cavail(n_c_types)
    real, parameter :: aerobic_max = 0.022 ! Maximum soil-moisture factor under ideal conditions



    enz=Chet*enzfrac

    !Available carbon is the amount that diffuses to enzyme site
    !Cavail=Ctotal*sol*D*theta**3
    !Available carbon is the amount that is not in saturated soil fraction
    !Cavail=Ctotal*(1.0-theta)
    !Power law dependence on soil moisture, see below

    Cavail=Ctotal
    IF(sum(Cavail).eq.0.0 .OR. theta.eq.0.0 .OR. enz.eq.0.0) THEN
        Resp=0.0
        return
    ENDIF

    Resp=Vmax(T)*theta**3*(Cavail)*enz/(sum(Cavail)*kC+enz)*max((air_filled_porosity)**gas_diffusion_exp,min_anaerobic_resp_factor*aerobic_max)

    !ox_avail=oxygen_concentration(Ox,sum(tempresp)/sum(Cavail)*theta*oxPerC)
    !print *,sum(tempresp)/sum(Cavail)
    !print *,ox_avail/(kO+ox_avail)
    !print *,tempresp*(1.0-theta)+theta*tempresp*ox_avail/(kO+ox_avail)

    !Assumes pores are either saturated or not oxygen limited, divided according to theta
    !Resp=tempresp*(1.0-theta)+theta*tempresp*ox_avail/(kO+ox_avail)
end function Resp


function Vmax(T)
    real,intent(in)::T
    real,parameter::Tref=293.15
    real,dimension(n_c_types)::alpha,Vmax

    alpha=vmaxref/exp(-Ea/(Rugas*Tref))
    Vmax=alpha*exp(-Ea/(Rugas*T))
end function Vmax


real function oxygen_concentration(Oxw,tw,kO,demand)
    real,intent(in)::Oxw,demand,tw,kO !Oxygen content of water, oxygen demand

    oxygen_concentration=0.5*(Oxw-tw*demand-kO + sqrt((tw*demand+kO-Oxw)**2 + 4*Oxw*kO))
end function oxygen_concentration

! =============================================================================
! prints on-line state of the carbon cohort
subroutine debug_cohort(c)
   type(litterCohort), intent(in) :: c
   __DEBUG___(c%litterC)
   __DEBUG___(c%protectedC)
   __DEBUG___(c%livingMicrobeC)
   __DEBUG___(c%CO2)
   __DEBUG___(c%originalLitterC)
   write(*,*)
end subroutine debug_cohort

! =============================================================================
subroutine debug_pool(p,tag)
   type(soil_carbon_pool), intent(in) :: p
   character(*), intent(in) :: tag

   integer :: i
   write(*,'(x,a)',advance='NO')trim(tag)
   __DEBUG5__(p%max_cohorts,p%n_cohorts,p%protection_rate,p%Qmax,p%dissolved_carbon)
!   do i = 1, p%n_cohorts
!      write(*,'(x,i3.3)',advance='NO')i
!      call debug_cohort(p%litterCohorts(i))
!   enddo
end subroutine debug_pool


subroutine print_cohort(cohort)
    type(litterCohort)::cohort

    WRITE (*,*) '----------------'
    WRITE (*,*) 'Original C =',cohort%originalLitterC
    WRITE (*,*) 'Unprotected C=',cohort%litterC
    WRITE (*,*) 'Living microbial C =',cohort%livingMicrobeC
    !WRITE (*,*) 'Microbial products =',cohort%microbProdC
    WRITE (*,*) 'Protected C =',cohort%protectedC

    !WRITE (*,*) 'Mineral complex C =',cohort%minC
    WRITE (*,*) 'CO2 =',cohort%CO2
    WRITE (*,*) 'Rtot =',cohort%Rtot
    !WRITE (*,*) 'Volume fraction =',cohort%cohortVolume
    WRITE (*,*) 'Sum of carbon =',cohortCsum(cohort)
    WRITE (*,*) '----------------'
end subroutine


!Check for carbon balance and invalid values
logical function check_cohort(cohort) result(cohortGood)
    type(litterCohort),intent(in)::cohort
    integer::n
    logical:: tempGood
    real :: cohortC
    ! ZMS
    real :: tol_roundoff = 1.e-11    ![kg C/m^2] tolerance for roundoff error in soil carbon numerics

    cohortGood=.NOT. ( &
    ! (min(cohort%originalLitterC,cohort%livingMicrobeC,cohort%CO2).lt.0) .OR. &
     (min(cohort%originalLitterC,cohort%livingMicrobeC,cohort%CO2).lt.-tol_roundoff) .OR. &
    is_nan(cohort%originalLitterC) .OR. &
    is_nan(cohort%livingMicrobeC) .OR. &
    is_nan(cohort%CO2)  &
    )

   cohortC=cohortCSum(cohort)
  !if(cohortC .GT. 0.0) cohortGood = cohortGood .AND. (abs(cohortC-cohort%originalLitterC)/cohortC) .LT. tol
  ! ZMS Make this an either relative or absolute balance check to allow roundoff error for infinitessimal pools in deep layers
    if(cohortC .GT. 0.0) cohortGood = cohortGood .AND. &
              (abs(cohortC-cohort%originalLitterC)/cohortC .LT. tol .OR. abs(cohortC-cohort%originalLitterC) .LT. tol_roundoff)

    DO n=1,n_c_types
        tempGood = .NOT. ( &
        (cohort%litterC(n).lt.0) .OR. &
        (is_nan(cohort%litterC(n))) .OR. &
        (is_nan(cohort%protectedC(n)))  &
        )
        cohortGood=cohortGood .AND. tempGood
    ENDDO

   ! IF(cohortBad) THEN
!        WRITE (*,*)'Cohort carbon pool error:'
!        call print_cohort(cohort)
!        WRITE (*,*)'Sum =',cohortCSum(cohort),'Difference =',cohortCSum(cohort)-(cohort%originalLitterC)

        !This should be changed to a model-friendly error
!#ifdef STANDALONE_SOIL_CARBON
!        !STOP
!#else
!        !call error_mesg('check_cohort','Error in cohort value',FATAL)
!#endif
!    !ENDIF
end function check_cohort



function cohortCSum(cohort,only_active)
    type(litterCohort),intent(in)::cohort
    logical,intent(in),optional::only_active
    real::tempSum
    real::cohortCSum

    logical::only_act

    only_act=.FALSE.
    if(present(only_active)) only_act=only_active


    tempSum=0.0
    tempSum=tempSum+cohort%livingMicrobeC
    tempSum=tempSum+sum(cohort%litterC)
    tempSum=tempSum+sum(cohort%protectedC)
    !tempSum=tempSum+sum(cohort%minC)
    if (.NOT. only_act) then
        tempSum=tempSum+cohort%CO2
    endif
    cohortCSum=tempSum
END FUNCTION


!Add a cohort to the soil carbon pool
subroutine add_cohort(pool,newCohort)
    type(soil_carbon_pool),intent(inout) :: pool
    type(litterCohort),intent(in) :: newCohort

    type(litterCohort),allocatable::newCohorts(:)

    if (.not.allocated(pool%litterCohorts)) then
       allocate(pool%litterCohorts(init_n_cohorts))
       pool%n_cohorts = 0
    endif

    if (pool%n_cohorts+1>size(pool%litterCohorts)) then
       ! not enough room for new cohort in the array of cohorts, reallocate
       ! cohort array: double its size
       allocate(newCohorts(2*size(pool%litterCohorts)))
       newCohorts(1:pool%n_cohorts) = pool%litterCohorts(1:pool%n_cohorts)
       call move_alloc(from=newCohorts,to=pool%litterCohorts)
    endif

    pool%litterCohorts(pool%n_cohorts+1) = newCohort
    pool%n_cohorts = pool%n_cohorts+1
end subroutine add_cohort


!Adds litter as a new cohort
subroutine add_litter(pool,newLitterC)
    type(soil_carbon_pool),intent(inout)::pool
    real,intent(in) :: newLitterC(n_c_types)

    type(litterCohort)::newCohort
    real::initialMicrobeC

    if(any(newLitterC < 0.0)) call error_mesg('add_litter','Some of new litter less than zero',FATAL)
    initialMicrobeC=sum(newLitterC*minMicrobeC)

    call initializeCohort(newCohort,litterInputC=newLitterC*(1.0-minMicrobeC),initialMicrobeC=initialMicrobeC)
    call add_cohort(pool,newCohort)

end subroutine add_litter

!Add carbon to all cohorts in pool, weighted by cohort size.
!For leaching or exudation that is spread evenly in soil instead of new cohort
subroutine add_carbon_to_cohorts(pool,litterC,protectedC,livingMicrobeC,CO2)
    type(soil_carbon_pool),intent(inout)::pool
    real,optional,dimension(n_c_types),intent(in)::litterC,protectedC
    real,optional,intent(in)::livingMicrobeC,CO2

    real,dimension(n_c_types)::litterCval,protectedCval
    real::livingMicrobeCval,CO2val
    real::totalCarbon,weight
    integer::ii

    livingMicrobeCval=0.0
    protectedCval=0.0
    litterCval=0.0
    CO2val=0.0

    IF (present(litterC)) litterCval=litterC
    IF (present(protectedC)) protectedCval=protectedC
    IF (present(livingMicrobeC)) livingMicrobeCval=livingMicrobeC
    IF (present(CO2)) CO2val=CO2

    IF (any(litterCval<0) .or. any(protectedCval<0) .or. livingMicrobeCval<0 .or. CO2val<0) &
          call error_mesg('add_carbon_to_cohorts','Carbon added less than zero',FATAL)

    call poolTotalCarbon(pool,totalCarbon=totalCarbon)
    if(.not.allocated(pool%litterCohorts)) call add_litter(pool,(/0.0,0.0,0.0/))
    DO ii=1,pool%n_cohorts
        !Deposited carbon is just weighted by total cohort carbon
        IF(totalCarbon.gt.0) THEN
           weight=cohortCsum(pool%litterCohorts(ii),only_active=.TRUE.)/totalCarbon
        ELSE
           weight=1.0
        ENDIF

        pool%litterCohorts(ii)%litterC=pool%litterCohorts(ii)%litterC+litterCval*weight
        pool%litterCohorts(ii)%protectedC=pool%litterCohorts(ii)%protectedC+protectedCval*weight
        pool%litterCohorts(ii)%livingMicrobeC=pool%litterCohorts(ii)%livingMicrobeC+livingMicrobeCval*weight
        pool%litterCohorts(ii)%CO2=pool%litterCohorts(ii)%CO2+CO2val*weight
        pool%litterCohorts(ii)%originalLitterC=pool%litterCohorts(ii)%originalLitterC+(sum(litterCval+protectedCval)+livingMicrobeCval+CO2val)*weight
    ENDDO

end subroutine


subroutine add_carbon_to_rhizosphere(pool,newCarbon,rhizosphere_frac)
    type(soil_carbon_pool),intent(inout)::pool
    real,intent(in)::newCarbon(n_c_types)
    real,intent(in)::rhizosphere_frac

    type(litterCohort)::rhizosphere
    real::litter_removed(n_c_types),protected_removed(n_c_types),liveMicrobe_removed
    integer::n

    call remove_carbon_fraction_from_pool(pool,rhizosphere_frac,litter_removed,protected_removed,liveMicrobe_removed)

    call initializeCohort(rhizosphere,litter_removed+newCarbon,liveMicrobe_removed)

    call add_cohort(pool,rhizosphere)
    do n=1,pool%n_cohorts
       IF (.NOT. check_cohort(pool%litterCohorts(n))) THEN
            WRITE (*,*) 'add_carbon_to_rhizosphere: Cohort',n,'of',pool%n_cohorts,'bad'
            call print_cohort(pool%litterCohorts(n))
            call error_mesg('add_carbon_to_rhizosphere','Bad cohort',FATAL)
       ENDIF
    enddo
end subroutine



!Remove a fraction of carbon from all cohorts.  For processes like leaching and bioturbation
subroutine remove_carbon_fraction_from_pool(pool,fraction,litter_removed,protected_removed,liveMicrobe_removed,litterMobility,protectedMobility,livingMicrobeMobility)
    type(soil_carbon_pool),intent(inout) :: pool
    real,intent(in) :: fraction
    real,intent(out) :: litter_removed(n_c_types),protected_removed(n_c_types),liveMicrobe_removed
    real,intent(in), optional :: litterMobility(n_c_types),protectedMobility,livingMicrobeMobility

    integer::ii
    real::temp,temp2(n_c_types)
    real::litterFactor(n_c_types),protectedFactor,liveMicrobeFactor

    litterFactor=1.0
    protectedFactor=1.0
    liveMicrobeFactor=1.0

    IF (present(litterMobility)) litterFactor=litterMobility
    IF (present(protectedMobility)) protectedFactor=protectedMobility
    IF (present(livingMicrobeMobility)) liveMicrobeFactor=livingMicrobeMobility

    if(fraction>1.0) call error_mesg('remove_carbon_fraction_from_pool','fraction > 1.0',FATAL)

    litter_removed=0.0
    protected_removed=0.0
    liveMicrobe_removed=0.0

    DO ii=1,pool%n_cohorts
        temp2=pool%litterCohorts(ii)%litterC*min(1.0,fraction)*litterFactor
        litter_removed=litter_removed+temp2
        pool%litterCohorts(ii)%litterC=pool%litterCohorts(ii)%litterC-temp2
        pool%litterCohorts(ii)%originalLitterC=pool%litterCohorts(ii)%originalLitterC-sum(temp2)

        temp2=pool%litterCohorts(ii)%protectedC*min(1.0,fraction)*protectedFactor
        protected_removed=protected_removed+temp2
        pool%litterCohorts(ii)%protectedC=pool%litterCohorts(ii)%protectedC-temp2
        pool%litterCohorts(ii)%originalLitterC=pool%litterCohorts(ii)%originalLitterC-sum(temp2)

        temp=pool%litterCohorts(ii)%livingMicrobeC*min(1.0,fraction)*liveMicrobeFactor
        liveMicrobe_removed=liveMicrobe_removed+temp
        pool%litterCohorts(ii)%livingMicrobeC=pool%litterCohorts(ii)%livingMicrobeC-temp
        pool%litterCohorts(ii)%originalLitterC=pool%litterCohorts(ii)%originalLitterC-temp
    ENDDO
end subroutine


subroutine remove_cohort(pool,num)
    type(soil_carbon_pool),intent(inout)::pool
    integer,intent(in)::num

    integer :: i
    !PRINT *,'Removing cohort',num
    IF (num > pool%n_cohorts) call error_mesg('remove_cohort','Tried to remove nonexistent soil cohort', FATAL)

    do i = num, pool%n_cohorts-1
       pool%litterCohorts(i) = pool%litterCohorts(i+1)
    enddo
    pool%n_cohorts = pool%n_cohorts - 1
    ! note that the cohort array does not shrink, because it is likely it'll
    ! need to be expanded in the future
end subroutine remove_cohort



subroutine combine_cohorts(cohort1,cohort2,result)
    type(litterCohort),intent(in)::cohort1,cohort2
    type(litterCohort),intent(out)::result

    real::totalVolume,Cbefore,Cafter

    call initializeCohort(result)

    result%litterC=(cohort1%litterC + cohort2%litterC)
    result%protectedC=cohort1%protectedC + cohort2%protectedC

    result%livingMicrobeC=(cohort1%livingMicrobeC + cohort2%livingMicrobeC)
    result%CO2=(cohort1%CO2 + cohort2%CO2)
    result%Rtot=(cohort1%Rtot + cohort2%Rtot)
    result%originalLitterC=(cohort1%originalLitterC + cohort2%originalLitterC)

    !Cbefore=cohortCsum(cohort1)*cohort1%cohortVolume+cohortCsum(cohort2)*cohort2%cohortVolume
    !Cafter=cohortCsum(result)*result%cohortVolume

    !IF (Cbefore.ne.Cafter) PRINT *,'Error in combining cohorts:',Cafter-Cbefore

end subroutine


!Returns the cohort with all fields multiplied by x
function multiply_cohort(cohort,x) result(result)
    type(litterCohort),intent(in) :: cohort
    real,intent(in) :: x
    type(litterCohort) :: result

    call initializeCohort(result,cohort%litterC,cohort%livingMicrobeC,cohort%protectedC,x,cohort%CO2)
    result%Rtot=cohort%Rtot*x

end function multiply_cohort


!Combine two soil_carbon_pool pools, with weighting
subroutine combine_pools(pool1,pool2,w1,w2)
    type(soil_carbon_pool),intent(in) :: pool1
    type(soil_carbon_pool),intent(inout) :: pool2
    real :: w1,w2,x1,x2
    integer::cc

    ! Make sure weights are normalized
    x1 = w1/(w1+w2)
    x2 = 1.0 - x1

    !First multiply existing cohorts by weighting
    DO cc=1,pool2%n_cohorts
        pool2%litterCohorts(cc)=multiply_cohort(pool2%litterCohorts(cc),x2)
    ENDDO

    !Then just add the cohorts in pool1 to pool2, with weights
    DO cc=1,pool1%n_cohorts
        call add_cohort(pool2,multiply_cohort(pool1%litterCohorts(cc),x1))
    ENDDO

    call cull_cohorts(pool2)
    pool2%dissolved_carbon=pool2%dissolved_carbon*x2 + pool1%dissolved_carbon*x1
end subroutine combine_pools


!transfers a fraction of one pool into another, preserving cohorts
subroutine transfer_pool_fraction(source, destination, fraction)
    type(soil_carbon_pool), intent(inout) :: source, destination
    real, intent(in) :: fraction

    integer::nn
    type(litterCohort)::transfer_cohort

    if(fraction>1.0) call error_mesg('transfer_pool_fraction','fraction > 1.0',FATAL)
    if(fraction<0.0) call error_mesg('transfer_pool_fraction','fraction < 0.0',FATAL)

    do nn=1,source%n_cohorts
        transfer_cohort=multiply_cohort(source%litterCohorts(nn),fraction)
        call add_cohort(destination,transfer_cohort)
        source%litterCohorts(nn)=multiply_cohort(source%litterCohorts(nn),1.0-fraction)
    enddo

    destination%dissolved_carbon=destination%dissolved_carbon+source%dissolved_carbon*fraction
    source%dissolved_carbon=source%dissolved_carbon*(1.0-fraction)

end subroutine

!Sums all cohorts in layer into a single cohort
type(litterCohort) function totalCarbonCohort(pool)
    type(litterCohort)::tempCohort
    type(soil_carbon_pool),intent(in)::pool
    integer::ii


    call initializeCohort(totalCarbonCohort)

    do ii=1,pool%n_cohorts
        call combine_cohorts(pool%litterCohorts(ii),totalCarbonCohort,tempCohort)
        totalCarbonCohort=tempCohort
    enddo
end function



subroutine poolTotalCarbon(pool,fastC,slowC,deadMicrobeC,liveMicrobeC,protectedC,&
        fast_dissolved,slow_dissolved,deadmic_dissolved,ncohorts,totalCarbon,&
        fast_protectedC,slow_protectedC,deadmic_protectedC)
    real, intent(out),optional::fastC,slowC,deadMicrobeC,liveMicrobeC,protectedC,totalCarbon,&
                                fast_dissolved,slow_dissolved,deadmic_dissolved,&
                                fast_protectedC,slow_protectedC,deadmic_protectedC
    integer, intent(out),optional::ncohorts
    type(soil_carbon_pool),intent(in)::pool

    type(litterCohort)::totalCohort
    integer::n

    totalCohort=totalCarbonCohort(pool)

    if (present(fastC)) fastC=totalCohort%litterC(1)
    if (present(slowC)) slowC=totalCohort%litterC(2)
    if (present(deadMicrobeC)) deadMicrobeC=totalCohort%litterC(3)
    if (present(liveMicrobeC)) liveMicrobeC=totalCohort%livingMicrobeC
    if (present(protectedC)) protectedC=sum(totalCohort%protectedC)
    if (present(fast_protectedC)) fast_protectedC=totalCohort%protectedC(1)
    if (present(slow_protectedC)) slow_protectedC=totalCohort%protectedC(2)
    if (present(deadmic_protectedC)) deadmic_protectedC=totalCohort%protectedC(3)
    if (present(fast_dissolved)) fast_dissolved=pool%dissolved_carbon(1)
    if (present(slow_dissolved)) slow_dissolved=pool%dissolved_carbon(2)
    if (present(deadmic_dissolved)) deadmic_dissolved=pool%dissolved_carbon(3)
    if (present(totalCarbon)) totalCarbon=sum(totalCohort%litterC)+sum(totalCohort%protectedC)+&
            totalCohort%livingMicrobeC+sum(pool%dissolved_carbon)

    if (present(ncohorts)) then
        ncohorts=0
        do n=1,pool%n_cohorts
            if (cohortCsum(pool%litterCohorts(n)).gt.0.0) ncohorts=ncohorts+1
        enddo
    endif
end subroutine poolTotalCarbon


function pool_total_carbon(pool) result(s)
   type(soil_carbon_pool), intent(in) :: pool
   real :: s

   integer :: i

   s = sum(pool%dissolved_carbon)
   do i = 1, pool%n_cohorts
      s = s + sum(pool%litterCohorts(i)%litterC(:))    &
            + sum(pool%litterCohorts(i)%protectedC(:)) &
            + pool%litterCohorts(i)%livingMicrobeC
   enddo
end function pool_total_carbon

!Combine cohorts in a pool, making sure the total number is less than the pool max cohorts
subroutine cull_cohorts(pool)
    type(soil_carbon_pool),intent(inout)::pool
    integer::n,ncombined,m,c1,c2
    real::Cbefore,Cafter,minDistance,r
    type(litterCohort)::tempCohort

    ncombined=0

    DO WHILE(pool%n_cohorts>pool%max_cohorts)
        !Find closest pair of cohorts
        minDistance = HUGE(1.0)
        do n = 1, pool%n_cohorts
           do m = n+1, pool%n_cohorts
              r = cohortDistance(pool%litterCohorts(n),pool%litterCohorts(m))
              if (r<minDistance) then
                 minDistance = r; c1=n; c2=m
              endif
           enddo
        enddo
        call combine_cohorts(pool%litterCohorts(c1),pool%litterCohorts(c2),tempCohort)
        !PRINT '("At cohort ceiling (",I3,"): combined cohorts ",I3," and ",I3,",difference=",F8.3)',pool%max_Cohorts,c1,c2,minDistance
        pool%litterCohorts(c1)=tempCohort
        call remove_cohort(pool,c2)
        ncombined=ncombined+1
    ENDDO

    Cafter=cohortCsum(totalCarbonCohort(pool))
    !IF(ncombined.gt.0) WRITE (*,*),'Combined',ncombined,'cohorts'
    !totalCombineError=totalCombineError+(Cafter-Cbefore)

end subroutine cull_cohorts


!Gives a logical answer of whether two cohorts can be combined, based on their similarity
logical function canCombine(cohort1,cohort2)
    type(litterCohort)::cohort1,cohort2
    real::sum1,sum2

    sum1=cohortCsum(cohort1)-cohort1%CO2
    sum2=cohortCsum(cohort2)-cohort2%CO2

    if(sum1.eq.0.0 .or. sum2.eq.0.0) then
        canCombine = .TRUE.
    else

    canCombine = (sum1.eq.0.0 .OR. sum2.eq.0.0) .OR. &
                (     Cpoolcomp(cohort1%litterC(1),sum1,cohort2%litterC(1),sum2,.FALSE.).LE.combineThreshold &
                .AND. Cpoolcomp(cohort1%litterC(2),sum1,cohort2%litterC(2),sum2,.FALSE.).LE.combineThreshold &
                .AND. Cpoolcomp(cohort1%litterC(3),sum1,cohort2%litterC(3),sum2,.FALSE.).LE.combineThreshold &
                .AND. Cpoolcomp(sum(cohort1%protectedC),sum1,sum(cohort2%protectedC),sum2,.FALSE.).LE.combineThreshold &
                .AND. Cpoolcomp(cohort1%livingMicrobeC,sum1,cohort2%livingMicrobeC,sum2,.TRUE.).LE.combineThreshold &
                )
     endif

end function


!Calculate a numerical estimate of the difference between two cohorts,
!for deciding which are the closest and can be combined
real function cohortDistance(cohort1,cohort2)
    type(litterCohort)::cohort1,cohort2
    real::sum1,sum2

    !sum1=cohortCsum(cohort1,.TRUE.)
    !sum2=cohortCsum(cohort2,.TRUE.)
    sum1=sum(cohort1%litterC)
    sum2=sum(cohort2%litterC)

    IF (sum1.eq.0.0 .or. sum2.eq.0.0) THEN
        cohortDistance=0.0
    ELSE

        cohortDistance=sqrt(Cpoolcomp(cohort1%litterC(1),sum1,cohort2%litterC(1),sum2,.FALSE.)**2 &
                    + Cpoolcomp(cohort1%litterC(2),sum1,cohort2%litterC(2),sum2,.FALSE.)**2 &
                    + Cpoolcomp(cohort1%litterC(3),sum1,cohort2%litterC(3),sum2,.FALSE.)**2 &
                    !+ Cpoolcomp(sum(cohort1%protectedC),sum1,sum(cohort2%protectedC),sum2,.TRUE.)**2 &
                    !+ (Cpoolcomp(cohort1%livingMicrobeC,sum1,cohort2%livingMicrobeC,sum2,.FALSE.))**2 &
                    )

    ENDIF

end function


!Calculate the difference between two cohort carbon pools.  Used to decide if the cohorts can be combined
real function Cpoolcomp(pool1,sum1,pool2,sum2,norm) result(compval)
    real, intent(in)::pool1,sum1,pool2,sum2
    logical,intent(in)::norm
    real::poolmean,pooldiff
    real::normval


    poolmean=0.5*(pool1/sum1+pool2/sum2)
    pooldiff=abs(pool1/sum1-pool2/sum2)

    IF (norm) THEN
        normval=poolmean
    ELSE
        normval=1.0
    ENDIF

    IF (poolmean.EQ.0.0) THEN
    compval=0.0
    ELSE
    compval=pooldiff/normval
    ENDIF
end function




subroutine tracer_advection(tracer_mass,flow,div,dz,wl,del_tracer,divergence_loss)
    real,intent(inout),dimension(:):: tracer_mass  !Per layer (not per unit water)
    real,intent(in),dimension(:)   :: flow  !Total flow, not flow rate [mm]
    real,intent(in),dimension(:)   :: div   !Horizontal divergence (layer total, not rate) [mm]
    real,intent(in),dimension(:)   :: dz    !Layer thickness
    real,intent(in),dimension(:)   :: wl    !water content [kg/m^2] by layer before Richards (1:num_l)
    real,intent(out),dimension(:)  :: del_tracer,divergence_loss !Change in tracer mass, and divergence part

    real,dimension(size(tracer_mass)) :: aaa,bbb,ccc,ddd    !Matrix coefficients for aaa*dx[i-1] + bbb*dx[i] + ccc*dx[i+1] = ddd
    real,dimension(size(tracer_mass)) :: u_minus,u_plus   !For weighting of flow upstream/downstream
    real,dimension(size(tracer_mass)) :: tracer_concentration ! [kg C/m^3 soil]
    integer::ll,nlayers
    real,dimension(size(tracer_mass)) :: flow_eff ! flow adjusted to be units of [m], weighted by 1/wl
    real, parameter                   :: minwl = 0.1 ![mm] minimum allowed wl

    nlayers=size(tracer_concentration)
    tracer_concentration=tracer_mass/dz
    ! ZMS
!    tracer_concentration(1) = tracer_mass(1)/dz(1)
!    do ll=2,nlayers
!       if (wl(ll-1) > 0.) then
!          tracer_concentration(ll) = tracer_mass(ll)/(wl(ll-1)/dens_h2o)
!       else
!          tracer_concentration(ll) = tracer_mass(ll)/dz(ll)
!       end if
!    end do

    ! Adjust flow to be in units of m, weighted by 1/wl
    ! wl of litter layer is undefined
    flow_eff(1) = flow(1)/dens_h2o ! This should be zero
    if (flow(2) > 0.) then
       flow_eff(2) = flow(2)/dens_h2o
       ! [m]       = [mm]/[mm/m]
    else
       flow_eff(2) = flow(2) * dz(2)/max(wl(1), minwl)
       ! [m]       = [mm]    * [m]  / [mm]
    end if
    do ll = 3, nlayers
       if (flow(ll) > 0.) then
          flow_eff(ll) = flow(ll) * dz(ll-1)/max(wl(ll-2), minwl)
       else
          flow_eff(ll) = flow(ll) * dz(ll)/max(wl(ll-1), minwl)
       end if
    end do

    u_minus = 1.
    where (flow.lt.0.) u_minus = 0.
    do ll = 1, nlayers-1
        u_plus(ll) = 1. - u_minus(ll+1)
    enddo

    !Top layer, uses upper bound concentration
    ! flow into top layer is zero
    ll=1
    aaa(ll)= 0.0 !flow(ll)*u_minus(ll)
    bbb(ll)= - flow_eff(ll+1)*(1-u_plus(ll)) - dz(ll)
    divergence_loss(ll)=max(div(ll),0.0)/dz(ll)/dens_h2o*tracer_concentration(ll)
    ! [kg/m^3]         =   [mm]         / [m]  / [mm/m] * [kg/m^3]
    ccc(ll)= -flow_eff(ll+1)*u_plus(ll)
    ddd(ll)= - tracer_concentration(ll)*(bbb(ll)+dz(ll)) - tracer_concentration(ll+1)*ccc(ll)
    do ll=2,nlayers-1
        aaa(ll)=flow_eff(ll)*u_minus(ll)
        bbb(ll)=flow_eff(ll)*(1-u_minus(ll)) - flow_eff(ll+1)*(1-u_plus(ll)) - dz(ll)
        divergence_loss(ll)=max(div(ll),0.0)/max(wl(ll-1),minwl)*tracer_concentration(ll)
        ! [kg/m^3]         =     [mm]       / [mm]              * [kg/m^3]
        ccc(ll)=-flow_eff(ll+1)*u_plus(ll)
        ddd(ll)=-tracer_concentration(ll-1)*aaa(ll) - tracer_concentration(ll)*(bbb(ll)+dz(ll)) - tracer_concentration(ll+1)*ccc(ll)
    enddo
    !bottom layer, flow out is zero
    ll=nlayers
    aaa(ll)=flow_eff(ll)*u_minus(ll)
    bbb(ll)= flow_eff(ll)*(1-u_minus(ll)) - dz(ll)
    divergence_loss(ll)=max(div(ll),0.0)/max(wl(ll-1),minwl)*tracer_concentration(ll)
    ccc(ll)= 0.0
    ddd(ll)=-tracer_concentration(ll-1)*aaa(ll) - tracer_concentration(ll)*(bbb(ll)+dz(ll))

    !Solve the linear algebra problem
    if(nlayers.gt.1) then
        call tridiag(aaa,bbb,ccc,ddd,del_tracer)
    else
        del_tracer=0.0
    endif


    del_tracer=del_tracer*dz
    tracer_mass=tracer_mass+del_tracer
    divergence_loss=divergence_loss*dz
    where(divergence_loss>tracer_mass) divergence_loss=tracer_mass

end subroutine tracer_advection




!Newer version using tracer_advection
!subroutine carbon_leaching(layers,flow,div,dz,dt,del_DOC)
!    type(soil_carbon_pool),dimension(:),intent(inout)::layers
!    real,dimension(:),intent(in)::flow,div,dz !flow in units of m3/m2, downward is >0
!    real,intent(in)::dt                       !In seconds
!
!    real,intent(out),optional::del_DOC(:,:)
!
!    real::litterbefore(n_c_types,size(layers))
!    real,dimension(n_c_types,size(layers)) :: dissolved_litter
!    real,dimension(n_c_types,size(layers)) :: div_loss_litter
!    real::d_litter(n_c_types,size(layers))
!    integer::ll,ii,badcohort,n
!    ! ZMS: for debug
!    integer:: i,j,k,face
!
!
!
!    do ii=1,n_c_types
!    litterbefore(ii,:)=layers(:)%dissolved_carbon(ii)
!    if(any(litterbefore(ii,:)<0)) then
!        print *,'Carbon flavor',ii
!        print *,litterbefore(ii,:)
!        call error_mesg('carbon_leaching','Dissolved litter < 0 (before advection)',FATAL)
!    endif
!
!    call tracer_advection(layers(:)%dissolved_carbon(ii),flow(1:size(layers)),div*dt*0,dz,d_litter(ii,:),div_loss_litter(ii,:))
!
!    if(any(layers(:)%dissolved_carbon(ii)<0)) call error_mesg('carbon_leaching','Dissolved litter < 0 (after advection)',FATAL)
!    enddo
!
!    if(present(del_DOC)) del_DOC=d_litter
!
!
!end subroutine carbon_leaching


subroutine carbon_leaching_with_litter(soil,leaflitter,woodlitter,flow,litterflow, &
      div,dz,wl,dt,del_soil_DOC,del_leaflitter_DOC, &
      del_woodlitter_DOC,div_DOC_loss,tiled,div_hlsp_DOC,surf_DOC_loss)
    type(soil_carbon_pool),dimension(:),intent(inout)::soil
    type(soil_carbon_pool),intent(inout)::leaflitter,woodlitter
    real,dimension(:),intent(in)::flow,div,dz,wl !flow and wl in units of mm, downward is >0
    real,intent(in)  :: litterflow  !Flow and divergence for litter layer
    real,intent(in)::dt                       !In seconds
    logical, intent(in) :: tiled ! flag for tiled hydrology
    real, optional, intent(in) :: div_hlsp_DOC(:,:) ! dim(n_c_types, num_l) [kg C/m^2/s] net divergence
                                                    ! loss from tile calculated in hlsp_hydrology
    real, optional, intent(out) :: surf_DOC_loss(n_c_types) ! [kg C/m^2] loss from top layer to surface runoff

    real,intent(out),optional::del_soil_DOC(:,:),div_DOC_loss(:,:),del_leaflitter_DOC(:),del_woodlitter_DOC(:)

    real::DOC(n_c_types,size(soil)+1)
    real,dimension(n_c_types,size(soil)+1) :: div_loss
    real::d_DOC(n_c_types,size(soil)+1)
    real,dimension(size(soil)+1) :: flow_with_litter, div_with_litter, dz_with_litter
    integer::l,ii,badcohort,n
    real::litterThickness,leaflitterTotalC,woodlitterTotalC,DOCbefore(size(soil)+1),leaf_DOC_frac
    real :: surf_DOC_loss_loc(n_c_types)
    real, parameter :: minwl = 0.1 ! [mm]

    !For now, use a mininum litter thickness of 5 mm
    call poolTotalCarbon(leaflitter, totalCarbon=leaflitterTotalC)
    call poolTotalCarbon(woodlitter, totalCarbon=woodlitterTotalC)
    litterThickness=max((leaflitterTotalC+woodlitterTotalC)/litterDensity,5e-3)

    flow_with_litter(1)=0.0
    !flow_with_litter(2)=litterflow
    flow_with_litter(2:size(flow_with_litter))=flow(1:size(flow_with_litter)-1)
    flow_with_litter=flow_with_litter/1000 !xz change the div unit from mm to m
    div_with_litter(1)=0.0
    div_with_litter(2:size(flow_with_litter))=div(:)*dt ! div is in mm/s
    div_with_litter=div_with_litter/1000 !xz change the div unit from mm to m
    dz_with_litter(1)=litterThickness
    dz_with_litter(2:size(dz_with_litter)) = dz(:)

    !flow_with_litter(:)=0.0
    !flow_with_litter(2:)=1e-3

    surf_DOC_loss_loc(:) = 0.

    do ii=1,n_c_types
       DOC(ii,1)=leaflitter%dissolved_carbon(ii)+woodlitter%dissolved_carbon(ii)

       if(DOC(ii,1)>0) then
           leaf_DOC_frac=leaflitter%dissolved_carbon(ii)/DOC(ii,1)
       else
           leaf_DOC_frac=0.0
       endif

       DOC(ii,2:size(soil)+1)=soil(:)%dissolved_carbon(ii)
       if(any(DOC(ii,:)<0)) then
           print *,'Carbon flavor',ii
           print *,DOC(ii,:)
           call error_mesg('carbon_leaching_with_litter','Dissolved carbon < 0 (before advection)',FATAL)
       endif
       DOCbefore=DOC(ii,:)
       if (is_watch_point()) then
          __DEBUG1__(DOCbefore)
       endif
       call tracer_advection(DOC(ii,:),flow_with_litter(:),div_with_litter(:), &
             dz_with_litter,wl,d_DOC(ii,:),div_loss(ii,:))
       if (is_watch_point()) then
          __DEBUG1__(DOC(ii,:))
       endif
       if(abs(sum(DOC(ii,:))-sum(DOCbefore)).gt.1e-10) then
           print *,'Flavor:',ii
           print *,'Before:'
           print *,DOCbefore
           print *,'After:'
           print *,DOC(ii,:)
           print *,'Difference:'
           print *,d_DOC(ii,:)
           print *,'Flow:'
           print *,flow_with_litter
           print *,'Total difference:',sum(DOC(ii,:))-sum(DOCbefore)
           call error_mesg('carbon_leaching_with_litter','Dissolved carbon not conserved',FATAL)
       endif
       if(any(DOC(ii,:)<0)) call error_mesg('carbon_leaching_with_litter','Dissolved carbon < 0 (after advection)',FATAL)

       if (tiled) then ! reset div_loss(ii,2:num_l+1) according to values calculated in hlsp_hydrology
          div_loss(ii,2:size(soil)+1) = div_hlsp_DOC(ii,:)*dt
          if (flow(1) < 0. .and. wl(1) > minwl) then  ! Add loss from top layer to runoff
             surf_DOC_loss_loc(ii) = -DOC(ii, 2) * flow(1) / wl(1)
             surf_DOC_loss_loc(ii) = min(surf_DOC_loss_loc(ii), DOC(ii,2))
          end if
          div_loss(ii,2) = min(div_loss(ii,2), DOC(ii,2) - surf_DOC_loss_loc(ii))
          do l=3,size(soil)+1
             div_loss(ii,l) = min(div_loss(ii,l), DOC(ii,l))
          end do
          ! Note: if these limits are imposed, there will be an imbalance between inter-tile fluxes
          ! that will be effectively rectified by subtracting from the flux to stream. In rare
          ! situations, that could lead to a negative stream DOC flux.
       end if
       DOC(ii,:)=DOC(ii,:)-div_loss(ii,:)
       DOC(ii,2)=DOC(ii,2)-surf_DOC_loss_loc(ii)

       leaflitter%dissolved_carbon(ii)=DOC(ii,1)*leaf_DOC_frac
       woodlitter%dissolved_carbon(ii)=DOC(ii,1)*(1.0-leaf_DOC_frac)
       soil(:)%dissolved_carbon(ii)=DOC(ii,2:size(soil)+1)

       if(present(del_soil_DOC)) del_soil_DOC(ii,:)=d_DOC(ii,2:size(soil)+1)
       if(present(del_leaflitter_DOC)) del_leaflitter_DOC(ii)=d_DOC(ii,1)*leaf_DOC_frac
       if(present(del_woodlitter_DOC)) del_woodlitter_DOC(ii)=d_DOC(ii,1)*(1.0-leaf_DOC_frac)

       if(present(div_DOC_loss)) div_DOC_loss(ii,:)=div_loss(ii,1:size(div_DOC_loss,2))
       if(present(surf_DOC_loss)) surf_DOC_loss(ii) = surf_DOC_loss_loc(ii)
    enddo


end subroutine carbon_leaching_with_litter


subroutine retrieve_DOC(soil, DOC, num_l)

   type(soil_carbon_pool),dimension(:),intent(in) :: soil ! soil carbon pointer
   integer, intent(in)  :: num_l ! number of soil layers
   real, intent(out)    :: DOC(n_c_types, num_l) ! [kg C/m^2] dissolved organic carbon
   integer :: l

   do l=1,num_l
      DOC(1:n_c_types,l)=soil(l)%dissolved_carbon(1:n_c_types)
   end do

end subroutine retrieve_DOC


subroutine litter_leaching(litterLayer,topSoilLayer,flow,dt,DOC_loss)
    type(soil_carbon_pool),intent(inout)::litterLayer,topSoilLayer  !Moves C between these layers
    real,intent(in)::flow       !Vertical water flow from litter layer to first soil layer
    real,intent(out)::DOC_loss(:) !Change in litter layer C from leaching (<0)
    real,intent(in)::dt  !dt in seconds

    real::litterDensity

    DOC_loss=0.0

    IF(flow.gt.0.0) THEN
        DOC_loss=litterLayer%dissolved_carbon*flow
            if(any(DOC_loss<0)) call error_mesg('litter_leaching','Litter flow < 0',FATAL)
            where(DOC_loss>litterLayer%dissolved_carbon) DOC_loss=litterLayer%dissolved_carbon
            litterLayer%dissolved_carbon=litterLayer%dissolved_carbon-DOC_loss
            topSoilLayer%dissolved_carbon=topSoilLayer%dissolved_carbon+DOC_loss
    ENDIF

    IF(.NOT. allocated(litterLayer%litterCohorts)) return

    IF(.NOT. check_cohort(litterLayer%litterCohorts(1))) THEN
         WRITE (*,*) 'LITTER_LEACHING: Cohort bad'
            call print_cohort(litterLayer%litterCohorts(1))
            call error_mesg('litter_leaching','Bad cohort after leaching',FATAL)
    ENDIF
end subroutine


#ifndef STANDALONE_SOIL_CARBON
subroutine adjust_pool_ncohorts(pool)
    type(soil_carbon_pool),intent(inout) :: pool

    !Remove cohorts if size is too large
    if (pool%n_cohorts.gt.soilMaxCohorts) call cull_cohorts(pool)

    !Add empty cohorts until size is correct
    do while (pool%n_cohorts.lt.soilMaxCohorts)
       call add_litter(pool,(/0.0,0.0,0.0/))
    enddo
end subroutine

#else

subroutine error_mesg(routine,message,level)
    character(len=*),intent(in)::routine,message
    integer::level

    if(level.eq.FATAL) then
    print *,'Error in routine ',routine,': ',message
    stop
    else
    print *,'Note from routine ',routine,': ',message
    endif
end subroutine

!Copied from soil_numerics (by sergey malyshev)
! ============================================================================
! given values of the triadiagonal matrix coefficients, computes a solution
subroutine tridiag(a,b,c,r,u)
  real, intent(in)  :: a(:),b(:),c(:),r(:)
  real, intent(out) :: u(:)

  integer :: j
  real :: bet, gam(size(a))

  ! check that the sizes are the same
  if(size(a)/=size(b).or.size(a)/=size(c).or.size(a)/=size(r)) &
       call error_mesg('tridiag','sizes of input arrays are not equal',FATAL)
  if(size(u)<size(a)) &
       call error_mesg('tridiag','size of the result is insufficient',FATAL)
  ! check that a(1)==0 and c(N)==0
  if(a(1)/=0.or.c(size(a))/=0) &
       call error_mesg('tridiag','a(1) and c(N) must be equal to 0',FATAL)
  ! decomposition and forward substitution
  bet = b(1)
  u(1) = r(1)/bet
  do j = 2,size(a)
     gam(j) = c(j-1)/bet
     bet = b(j)-a(j)*gam(j)
     if(bet==0) &
          call error_mesg('tridiag','system is ill-defined',FATAL)
     u(j) = (r(j)-a(j)*u(j-1))/bet
  enddo
  ! backward substitution
  do j = size(a)-1,1,-1
     u(j) = u(j)-gam(j+1)*u(j+1)
  enddo
end subroutine tridiag


#endif

! pgi: does not have bult-in isNaN function
! gfortran: versions <5 do not appear to support ieee_arithmetic module

! note that there is a danger that the compilr optimizes the comparison away
! but I don't see how
logical elemental function is_nan(x)
   real, intent(in) :: x
   is_nan = (x/=x)
end function is_nan

end module soil_carbon_mod
